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ABSTRACT 



Beginning with a relativistic action principle for the irrotational flow of collisionless mat- 
ter, we compute higher order corrections to the Zel'dovich approximation by deriving a 
nonlinear Hamilton-Jacobi equation for the velocity potential. It is shown that the velocity 
of the field may always be derived from a potential which however may be a multi-valued 
function of the space-time coordinates. In the Newtonian limit, the results are nonlocal 
because one must solve the Newton- Poisson equation. By considering the Hamilton-Jacobi 
equation for general relativity, we set up gauge-invariant equations which respect causal- 
ity. A spatial gradient expansion leads to simple and useful results which are local — they 
require only derivatives of the initial gravitational potential. 
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1. INTRODUCTION 



Recent redshift surveys of galaxies suggest that sheet-like structures are quite com- 
mon in our Universe. Within a distance of 100/i _1 Mpc, a Great Wall appears in the Center 
for Astrophysics slice of the Universe, De Lapparent, Geller & Huchra (1986). A pencil 
beam survey, Broadhurst et al (1990) is consistent with the hypothesis that similar sheets 
exist out to a redshift of z ~ |. The Zel'dovich (1970) approximation plays a fundamental 
role in describing the formation of sheets that arise from the nonlinear evolution of irrota- 
tional collisionless dust (for a review, consult e.g., Shandarin & Zel'dovich (1989). In an 
attempt to better understand the formation of cosmic structure, we show how to compute 
higher order terms in the Zel'dovich approximation. We employ two principal tools: (1) 
a relativistic action principle for dust interacting with gravity, and (2) Hamilton-Jacobi 
(HJ) theory used in conjunction with a Taylor series expansion. 

Dust with rest mass m which flows irrotationally may be described using an action 
principle for a dust field %, 

Idust = J d 4 xv^{-±p{l + g^ X ^X,v/m 2 ) }. (1.1) 

The four-velocity 

U» = -g' u 'x,u/m. (1-2) 

is proportional to the gradient of the potential x- The energy density p is a Lagrange 
multiplier which ensures that the four- velocity squared is —1: 

£W M = -1. (1.3) 
Variation with respect to x leads to the equation of conservation 

(pr% = o, 

whereas variation with respect to the 4-metric g^ u leads to the usual stress-energy tensor: 

This action principle (1.1) can be used to clarify the analysis of cosmo logical systems. It 
is relativistic and nonlinear. It will be the focal point of our discussion. 

The action for gravity is the standard Einstein-Hilbert form 

"-^gravity ~ 

where is the Ricci scalar of the 4-metric. Here we have taken 8ttG = 1. We shall 
also assume that m = 1 which can always be arranged by rescaling: x ~ ¥ m X- The action 
principle for dust is a special case of a general class of fluids considered by Schutz (1970, 
1971); see also Salopek & Stewart (1992, 1993) who obtain exact solutions for a two-fluid 
system of black-body radiation and dust. 
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Hamilton-Jacobi theory facilitates the generation of solutions for the corresponding 
field equations. In section 2, we write down the non-relativistic HJ equation for a cos- 
mological dust field, Kofman 1991). The generating function S is essentially a nonlinear 
velocity potential which is closely related to what is directly observable in cosmology (see 
e.g., Dekel et al (1990)). (We assume throughout that, as seen in a Lagrangian descrip- 
tion, the velocity is derivable from a potential. However when caustics form S may be 
a multivalued function of the space-time coordinates. In an Eulerian description several 
irrotational flows are superposed and the mean, Eulerian, flow is no longer irrotational.) 
The Zel'dovich approximation is the first two terms in our Taylor series expansion. Our 
technique clarifies the method employed by Buchert (1989, 1992, 1993) and Moutarde et 
al (1991). However, the higher order terms in the series are difficult to treat analytically 
because one must solve the Newton-Poisson equation, which is non-local. 

In section 3, we bypass some of the difficulties introduced by the Newton-Poisson 
equation by employing the relativistic Hamilton-Jacobi equation which is a functional dif- 
ferential equation. Although the formation of structure in a cosmological setting is basi- 
cally non-relativistic, the relativistic equations are useful because they maintain causality 
(which is not true for the Newtonian theory) . The generating functional S is solved using 
a spatial gradient expansion. By choosing our time hypersurface to be a slice where the 
dust-field is uniform, we are able to integrate the resulting equations using an iterative 
technique. In this way, we solve the energy and the momentum constraints of general rel- 
ativity. The energy constraint is essentially the G® Einstein equation; it is the relativistic 
generalization of the Newton-Poisson equation. The G® Einstein equation yields the mo- 
mentum constraint which ensures that the theory is invariant under reparameterizations 
of the spatial coordinates ( "gauge-invariance" ) . 

Our spatial gradient expansion is similar to solving Einstein's equations using a Tay- 
lor series, an approach that was initially pioneered by Lifshitz & Khalatnikov (1964). 
They were interested in describing the initial singularity of the big bang (see also Landau 
& Lifshitz (1975)). However their program was incomplete in that they did not explicitly 
solve the momentum constraint equation, which we solved recently under quite general 
conditions (Croudace, Parry, Salopek & Stewart (1994) hereafter known as CPSS; see also 
Salopek, Stewart, Croudace & Parry (1994)). Using the first two terms in our improved 
gradient expansion, we generalized the Zel'dovich approximation to a relativistic frame- 
work. We also investigated higher order terms which involve only derivatives of the initial 
gravitational potential — the results are local in nature. In this way, we hope to better 
understand the formation of sheet-like structures in the Universe. 

Matarrese et al (1993) considered an alternative relativistic approach to the formation 
of cosmic structure. They simplified and approximated Einstein's equations by neglecting 
the magnetic part of the Weyl tensor. The resulting equations were local in that they did 
not require any information about neighboring points. They determined an exact solution 
and they also integrated the equations numerically. CPSS demonstrated that the exact 
solution of Matarrese et al (1993) corresponded to the exact Szekeres solution of Einstein's 
equations. They showed further that this solution was unstable. Bertschinger & Jain 
(1994) suggested that the outcome of this instability was cigar shaped objects which were 
more generic than sheets. Counter to one's intuition, they also claimed that nonlinear 
structures can form in void regions. However, although this approach looks promising, 
it is still not clear when the approximation of neglecting the magnetic part of the Weyl 
tensor breaks down. 
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In section 4, we compare the results of our gradient expansion with two exact solutions 
of Einstein's equations, namely planar geometries and spherical geometries (Tolman-Bondi 
solutions). In general, the higher order corrections to the Zel'dovich are significant. In 
section 5 we compare our results with the second order theory of Buchert & Ehlers (1993), 
and in the final section we discuss our results and state our conclusions. 



2. EVOLUTION OF NONLINEAR VELOCITY POTENTIAL 



The Newtonian limit is essentially a non-relativistic approximation to general rel- 
ativity. In special relativity the energy of a particle with rest mass m is related to its 
momentum through 

E = m\J\ +p 2 /m 2 . 

The non-relativistic limit is the first two terms in the Taylor series (in p/m) which yields 

2 

In this section, we will use a similar technique to derive the cosmological HJ equation for 
collisionless dust. 



2.1 Non-relativistic Hamilton- Jacobi Equation 



In a cosmological setting, the dynamics of the matter is contained in the scalar 
evolution equation (1.3) for the dust field: 

(x - N* X J/N = sjl + x^ 1 - (2.1a) 

Here we have decomposed the metric using the Arnowitt-Deser-Misner (ADM) formalism 
(see, e.g., Misner et al (1973)): 

g 00 = -N 2 +Y J N,N J , g 0i = g i0 = N u g l3 = 7y . (2.16) 

N, Ni and 7^ are the lapse function, the shift function, and the 3-metric, respectively. A 
vertical bar | denotes the covariant derivative with respect to the 3-metric 7^. 

We now assume that spatial gradients are small and we approximate the right-hand 
side of eq.(2.1a) by the first two terms of a binomial series: 1 + §X|iX'*- We also work with 
a longitudinal gauge metric: 

N = l-$ H (t,x), N l = 0, 7ij =a 2 (t)[l + 2$ H (t,x)}5 lj; (2.2) 

here a(t) = t 2 / 3 is the scale factor, and $#(£, x) is the gauge-invariant variable of Bardeen 
(1980) (see also Stewart (1990)). The x l are comoving spatial coordinates. In cosmological 
problems which are described using longitudinal gauge, the deviations of the metric from 
Friedman- Robertson- Walker are of the order of 10~ 5 , and a linear treatment of the metric 
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is justified (this is not true in alternative gauges such as synchronous gauge; see section 
3). Next, we define the generating function S through 



X (t,x) = t-S(t,x) (2.3) 

which leads to the non-relativistic Hamilton- J acobi equation for a particle moving in a 
time-dependent potential — <&#(£, x): 

0=— + — — —-$#. (2.4a) 
dt 2a 2 dx l dx l v ; 

All the terms that were dropped are negligible. The nonrelativistic HJ equation was 
derived earlier by Kofman (1991) who noted that it was similar to the eikonal equation in 
geometric optics. If for the sake of illustration the potential is fixed, then this equation 
will exhibit caustics — i.e., S will be multivalued. However, this does not prevent one from 
solving the equation (Goldstein (1981)). 

The definition of the velocity potential and the four-velocity, 

U - ^ - 9 X ,u, 

(where r is proper time) implies that the coordinate velocity of a particle is related to S 
through 

dx l 1 dS 

-* =;?«?• (2 ' 4 ") 

Hence, the generating function S is essentially a velocity potential. Finally, the standard 
analysis of small perturbations of the metric leads to the Newton-Poisson relation 

-Av^ H =(p-A). (2.4c) 

The right-hand side is just the density perturbation 5p co on a comoving slice. In general, 
S and hence p will be multivalued functions, in which case the right-hand is summed over 
all particles with the same spatial coordinate x 1 . Finally, far within the Hubble radius, it 
is sufficient to approximate x ~ t in the last equation. 

It is easy to see that these equations are equivalent to the standard Newtonian equa- 
tions. If we differentiate eq.(2.4a) with respect to , one finds that 

d (dS\ oV d*S d$H . , 



dt \dxi J dt dx l dxi dxi 
which may be rewritten using a convective derivative, 



d _ d dx l d , , 

di = dt + ltt'dx~ i ' ( ' ' 
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r l = ax 1 



leading to 

° = dt (dxl) ~ ~dxT ' ^ 2 ' 7 ^ 
The physical displacement of a particle is denoted by 

and its velocity is 

d , d ,• dx 1 

— (ax) = -r + a—— . 
dt a dt 

In order to compute the peculiar velocity v l , one subtracts the velocity, (d/a)r\ of the 
Hubble flow, 

d , a a 1 dS 

v = — (ax) r = --— . 

dV ' a adx 1 

Substituting in eq.(2.7), one obtains the standard cosmological equation (Peebles (1980)): 

The density is computed from the equation of continuity which implies that 

(dx 1 \ 
g-j , (2-9) 

where qi denotes the initial (Lagrangian) coordinate, and po(t) = 4/(3t 2 ) is the energy 
density at early times. 

It is reasonable to assume that one does not commit any gross error in assuming that 
the velocity is derived from a potential S. However, one must be aware that the velocity 
potential may become multivalued, and even discontinuous. It is interesting to compare 
our equations with the standard linear perturbation analysis of the velocity potential where 
one neglects the second term in eq.(2.4a) which is a quadratic perturbation. This truncated 
system is a poor approximation for the formation of cosmological objects because the virial 
theorem implies that for a bound system the kinetic energy is of the order of the potential 
energy. 

Equations (2.4) are cumbersome to solve directly since we must solve S as a function 
of x % which in turn is a function of the Lagrangian variables qK We will now rewrite 
them in a form that is more in character with the HJ formalism that we are developing. 
We begin by inverting the spatial variables: we assume that the Lagrangian coordinates 
qi are functions of the Euler coordinates x l . The evolution equations for the Lagrangian 
coordinates are easy to write since they are constants in time: dq^ /dt = 0. As a result, 
one can write eqs.(2.4) in the form 

n dS 1 8S dS ^ 

° = ^7 + TT-o^TTr--$H, (2.10a) 



dt 2a 2 dx 1 dx 1 
dqi 1 dSdqi 

0= ^F + «^^' (2 - 106) 



2 4 ( n ^ (dq> 

a, 



^^Idetl^]-!). (2.10c) 
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In this way all the dependent variables are functions of the Eulerian coordinates x l . 



2.2 Series Solution of Non-relativistic HJ Equation 



We will now solve these cosmological equations by making a Taylor series expansion 
in the longitudinal time variable t: 



S(t, x) = 1 (x) + 1 5 / 3 (x) + t 7 ' 3 ( x ) + .... 
q j (t, x) = x j + t 2 ' 3 (x) + t 4/S qMi (x) + . . . , 
& H (t,x) = $%\x)+t 2/3 §%\x)+t A/3 <S>%\x) + ... . 



(2.11a) 
(2.116) 
(2.11c) 



Collecting common terms we find 



o = s<°) - <&g> , 








H ' 



5 g(1) | I dSMdSW 
3 2 cte* dx l 

%i , 



dx l 



(2.12a) 
(2.126) 

(2.12c) 



= 



V w +— , 

3 dxi ' 

4 (1)J a^w ggw ggWf 

+ ^7— + 7— + 



3 

2qW 



dxi 



dx l dx l 



dx l dx l 



(2.13a) 
(2.136) 
(2.13c) 



V 2 $g } = 



V 2 d>£ } = 



2dq 



(0)i 



3 cV 
2 cV ^ 



3 3 dx i dx? 



3 <9a>> cV 



(2.14a) 
(2.146) 



The solution for the zeroth order terms is 



s<°> = *g> = /(*) , <z (o)j = -?!^. 



(2.15a, 6) 



The arbitrary function of the spatial coordinates, —f(x), may be interpreted as the initial 
gravitational potential. The first order terms are expressed in terms of the non-local 
quantity W 



W 



d 2 f d 2 f 
dx k dx l dx k dx l 



- (v 2 /) ; 



(2.16a) 
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as 



SW(x) 



l^ j (x) 



3 df df 9 Jir 
- — -\ W 

4 dx i dx i 14 ' 

3 df df 15 
— H W 

4 dx i dx i 14 ' 



d 

dxi 



9 df_ df_ _ 2J_ w 

8 dx i dx i 56 



(2.166) 
(2.16c) 
(2.16d) 



The Zel'dovich approximation is obtained by solving for x l in favour of the Lagrangian 
coordinate q k using an iterative technique in eq.(2.11b): 



2 ag 1 56 dq 1 



d 2 f d 2 f 
QqkQqi QqkQqi 



(v 2 /) ; 



(2.17) 



Here / and its derivatives are evaluated at the point g J (which is independent of time t) 
rather than at x l which was the case for eqs.(2.16). The first two terms are the usual 
Zel'dovich approximation, and the third term is essentially the improvement suggested 
by Buchert and Ehlers (1993). Here HJ methods simplify the derivation. For a planar 
gravitational system, f(q l ) = f(q 3 ), where the initial gravitational potential is a function 
only of the third spatial coordinate, the matrix of second derivatives for / is diagonal with 
a single non-zero component 



d 2 f 
dq l dq J 



diag 



0,0, 



d 2 f 
dq 3 dq 3 



(2.18) 



and the last term in eq.(2.17) vanishes; hence the Zel'dovich approximation is essentially 
exact. The density is found using eq.(2.9). For a planar geometry, it is 



P = 



3t 2 \ l + ^t 2 / 3 (d 2 f/dq 3 dq s )\ ' 



(2.19) 



In general, further analytic progress is hampered by the non-local operator V 2 
appearing in eq.(2.17). 



3. RELATIVISTIC HAMILTON- J ACOBI THEORY 



Using a general relativistic formulation, we show how to compute the higher order 
terms in the Zel'dovich approximation which describes cosmological collapse. We utilize 
the Hamilton- Jacobi equation for general relativity which is a functional differential equa- 
tion that was first written by Peres (1962). We evolve the 3-metric in a spatial gradient 
expansion which is gauge-invariant. Our method is a significant advance over a Newtonian 
theory because it is local at each order. Using the improved Zel'dovich approximation, we 
compute the epoch of collapse. 
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Although the HJ formalism for general relativity can be applied to arbitrary time- 
hyper surfaces, in this section we will restrict ourselves to comoving surfaces where the 
velocity potential x is uniform. The line element is then 



ds 2 = -dx 2 + lij (x, q)dq l dq 3 . 

This choice is referred to as comoving, synchronous gauge, in which case the 3-metric 7^ 
has physical significance: it gives the physical distance, ds = ('jijdq l dq : >) 1 / 2 , between two 
comoving observers separated by an infinitesimal comoving distance dq l . x is related to 
the longitudinal time variable t through eq.(2.3). 

The relativistic Hamilton- J acobi equation for the restricted generating functional 
S[lij{x)\x] is 



dS f 7 o f _ 1/2 5S 5S 
= - — h / dTx { 7 1/2 



dx J I ^ij(x) S-f M (x) 



[2 ljk {x) lu {x) -7y(a07fci0*)] " l^ 1/2R ] ■ 



It is a functional differential equation. In addition one must satisfy the the momentum 
constraint 

/ 5S \ 5S , nn . 

= -2 7ife - — +-— _ 7IM , 3.2 

which states that the restricted generating functional «5[7y(x)|x] remains invariant under 
reparametrizations of the spatial coordinates. Eq.(3.1) is essentially the Gq Einstein equa- 
tion — it is the relativistic generalization of the Newton- Poisson equation (2.4c). Eq.(3.2) 
is the G® Einstein equation, which enforces gauge invariance. If eqs.(3.1) and (3.2) have 
been solved for the generating functional S, then all of Einstein's equations are satisfied 
provided that the 3-metric evolves according to 

^ = 2 7 " 1/2 ~r~~7~T ( 27*7i« - HHkl) ■ (3.3) 



dx $7ki(x) 
This is a first-order nonlinear differential equation for 7^ . 

3.1 Gradient Expansion for Generating Functional 

We look for a solution of the relativistic HJ equation (3.1) in the form of a spatial 
gradient expansion 

S = S(°) +5(2) +..., (3.4) 



where 



S<°> = - 2 J dW /2 H( X ) , (3.5a) 
S^=jdW /2 J(x)R, (3.56) 
S (4) = J d 3 £ 7 1/2 [L(X)R 2 + M( X )R ij Rij] ■ (3.5c) 



contains no spatial gradients whereas contains two spatial gradients, and so on. 
These functionals are invariant under transformations of the spatial coordinates since they 
are constructed using the invariant volume d 3 x , y 1 / 2 and the Ricci tensor R^i of the 3- 
metric 7^. H, J, L, M are arbitrary functions of the velocity potential x which are chosen 
so as to satisfy the relativistic HJ equation order by order: 



2 dH 

H 2 H — = , (zeroth order) (3.6a) 

3 dx 

^ + JH - ]- = , (second order) (3.66) 
ax 2 



(fourth order) (3.6c) 
(fourth order) (3.6d) 

As a result, the functional differential equation (3.1) has been reduced to a series of ordinary 
differential equations (3.6) which possess the trivial solution: 



dL 

dx 
dM 

d,\> 



LH 



:J = 0, 



- MH + 2 J 2 



2 

H = — , (zeroth order), (3.7a) 
3 

J=— X, (second order) , (3.76) 

L= 2^0 X3 ' M = ~Eo x3 ' ( fourthorder )- ( 3 - 7c ' rf ) 

The constant of integration in (3.7a) can be suppressed by choosing the origin of x a P _ 
propriately. The constants of integration for (3.7b-d) become irrelevant at large x an d 
have been ignored. Using HJ theory, we have thus been able to solve the constraints that 
appear within Einstein's equations. 



3.2 Gradient Expansion for 3-metric 



We solve the evolution equation for the 3-metric by using an iterative technique. If we 
neglect all second order spatial gradients and replace S by the zeroth order approximation 
iS(°) in eq.(3.3) we find that the 3-metric evolves according to 

7il ) (x,9)=X 4/3 %(?), (3.8) 

which is accurate to first order in spatial gradients. This is the long- wavelength approx- 
imation. The 'seed metric' kij(q) is an arbitrary function that is independent of time; 
it describes the initial fluctuations whose wavelengths are larger than the Hubble radius. 
This result is very old — it was known to Lifshitz & Khalatnikov (1964). 
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We thus derive the higher order corrections in a systematic way. For example, accu- 
rate to third order in spatial gradients, we have found that 



7gW) = X 4/3 



(3rd order) (3.9) 



and to fifth order, we obtained 

(5), , — 9 

81 



%\x,q) = x 4/3 k lJ + —r 



20' 



Rk, 



13 



AR l3 



+ 



(5th order) 



(3.10) 



350 



X 



8/3 



k 



AR m Ri m + —R z + qR-u ,k ) + ~~ lORRij + 17R l iRij 



89 



32 



8 



8 



5n » 

2 i r> k 



Here Rij(q) is the Ricci tensor of the seed metric, and a semi-colon (;) denotes its cor- 
responding covariant derivative. (For further details, see CPSS.) Using an alternative 
method, these expressions have been verified by Comer et al (1994). 



3.3 Improved Gradient Expansion for 3-metric 



Unfortunately, as time increases, the third order expression, eq.(3.9), leads to non- 
sensical results: the determinant of the 3-metric can actually become negative. A similar 
problem occurs when one naively approximates a non-negative function cos 2 (;r) ~ 1 — x 2 
with the first two terms of a Taylor expansion — the approximate function falls below 
zero when x > 1. One can easily remedy this problem by expanding cos(a;) ~ 1 — x 2 /2 , 
and then approximating cos 2 (a;) by the square of this result — this technique guarantees 
a positive result. In an analogous way, we can improve the expansions (3.9) and (3.10) by 
expressing the results as a 'square.' The improved 3-third order result 



i 3) (x,?) = x 4/3 \ k, + ^ ; 
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Rk u - ARu 



(3.11) 

yields the relativistic generalization of the Zel'dovich approximation (see section 3.4). The 
improved fifth order result is 



^\x,<l)=X 4/3 \ku + ^X 2/3 



40' 



Rku - ARu 



+ X 4/3 F tl k 



kjm + 4Q X 



2/3 



Rkj m ARj m 



+ X 4/3 F, 



jm ( 7 



(3.12a) 



where 



13 5600 



-32R lm R lm + 5R, k ' k + y£ 2 



-66^, +108^^ ;-20R l7 . k ' k +5R 



(3.126) 
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In comoving synchronous gauge, conservation of particle number implies that P7 1 / 2 is 
independent of time, and hence the number density p is given by 

p(x,4) = ^ 1/2 (4)7- 1/2 (x,9). (3-13) 

More details are given in Salopek & Stewart (1992), CPSS and Parry et al (1994). 

The argument for using the 'square' trick is heuristic, and it can placed on firmer 
theoretical grounds by showing that it leads to an exact solution of Einstein's equations. 
Our gradient expansion solution is essentially a Taylor series that is accurate for small \. A 
'priori there is no reason that it should work well in the regime where pancake structures are 
forming. However, after choosing an appropriate seed, the expressions for the 'square' agree 
precisely with the exact Szekeres (1975) solution which describes axisymmetric pancake 
formation (see CPSS). Hence, we have essentially fitted the late-time behavior to the exact 
solution. This method is similar to the application of WKB approximation in quantum 
mechanics. There, one matches two WKB solutions across the classical turn-around point 
by using the exact Airy solution for a linear potential. 



3.4 Modeling the Observable Universe 



Further simplifications arise when we wish to model our observable Universe. Al- 
though gravitational radiation arising from inflation (Salopek 1992) can indeed contribute 
significantly to the cosmic microwave background anisotropy observed by the COBE satel- 
lite, its interaction with cold-dark-matter is negligible during the matter-dominated era. 
For the purposes of describing the formation of galaxies, it is adequate to write the initial 
seed metric in the isotropic form: 

The arbitrary function f(q) is Bardeen's gauge-invariant variable at early times, and it 
can be interpreted as the negative of the initial gravitational potential. Whereas all of 
the expressions in section 3.3 were invariant under the group of spatial coordinate trans- 
formations, eq.(3.14) is invariant only under the subgroup of spatial rotations. We have 
thus sacrificed some of the gauge-invariance in favor of expressing the seed-metric in a very 
simple form involving the single function f(q). (However, in the derivation of eqs.(3.11), 
(3.12), it is very useful to assume that all quantities are fully gauge-invariant, otherwise 
the expressions become unwieldy.) One may readily compute the Ricci tensor for the above 
seed metric: 

RiM = 3 (-/,y - 8vf\i) + -q (/.,/•., " f,if'%) • (3-15) 

Our results may be simplified even further if we make one additional approximation: 
at each order in the gradient expansion, we will retain only the lowest order terms in / 
since / ~ 10~ 5 is much smaller than unity. For example, we obtain 

7^(x,<z) = X 4/3 % (first) (3.16) 
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at first order, 



% (3) (x,<z) = x 4/3 



6u + |x 2/3 /, 



, (improved 3rd) 



(3.17) 



at improved third order, and 



4/3 



Su + ^X 2/3 f,u + X 4/3 Fu 



4 + |x 2/3 /^ + x 4/3 ^- 



at improved fifth order where 

27 



(improved 5th) 
(3.18a) 

(3.186) 



Fa \ {f,uf> l ,j - t l ,if,ij) + s„ (if 1 .,) 2 - /•'"' /./,„); • 

It is important to note that the number density at the improved third order is exactly the 
result quoted by Zel'dovich: 



p (3) (x,g) 



3x 2 



/det 



(3.19) 



Hence we are justified in claiming that our improved third order spatial gradient approxi- 
mation reproduces the usual Zel'dovich approximation (CPSS, Salopek et al (1994)). 



3.5 Calculation of Collapse Epoch 



In order to understand the consequences of the improved fifth order results, we will 
compute the time when the density becomes infinite. 

For a given spatial point, we will rotate our coordinates so that the matrix of second 
derivatives is diagonal, 



f,ij — 



d 2 f 
dq i dq^ 



diag[-di,-d 2 ,-d3] , 



(3.20) 



Fortunately, the 3-metric at improved fifth order, eq.(3.18), is also diagonal, which gives 
it a huge advantage over the original form eq.(3.12). The density at improved fifth 
order becomes infinite when 



det 



8u + h 2: \r.u + 1\' [4 {f,uf' l ,j - f' l ,if,ij) + kj (if- 1 .,) 2 - f Jni fj,„). 



(3.21) 

If di is the largest positive eigenvalue of —f,ij, then this occurs when the scale factor 
(a = x 2 / 3 ) is given by 



a (5) - — / 
acoU ~ 3d, 1 



(3.22) 
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which is a function of q. In Fig. 1 we plot the collapse time for various values of the order 
pair (d,2/di, ds/di) which describe various ellipsoid configurations of the quadratic form 
f j ijAq t Aq : i = 1. (1,1) is a sphere, (1,0) and (0,1) are cylinders, and the line o^/di = d%/d\ 
represents a class of oblate spheroids. 

The origin (0,0) corresponds to the situation when d\ » \d2\,\ds\ which describes 
a sheet or plane. As expected from the discussion at the end of section 2.2, the standard 
Zel'dovich approximation (improved 3rd order) is essentially exact for this case. However, 
when all the eigenvalues are equal, ^1 = ^2 = ^3, describing a locally spherical system, 

then cicoll/ a< coll = 0-755 represents a 24% decrease over the improved third order result. 
The fifth order correction can be quite significant although it probably breaks down before 

one reaches ci^ n /af^ n = 2 shown in Fig. 1. 
4. COMPARISON WITH EXACT SOLUTIONS OF EINSTEIN'S EQUATIONS 



In CPSS, we compared our improved spatial gradient expansion with the Szekeres 
solution. We will now compare with two well-known exact solutions: planar and spherical 
geometries. 



4.1 Planar Geometries and the Zel'dovich Approximation 



For a planar gravitational system, the Zel'dovich approximation is essentially exact 
until caustic formation provided that the initial gravitational potential f(z) satisfies, 

f(z) « 1, (4.1a) 

and 

where z denotes the spatial coordinate. This result is well-known for the Newtonian theory 
(see, i.e., Kofman (1991)) but here we will use the methods of general relativity. For a 
perturbation f(z) = focos(kz), with fo < 10 -5 consistent with observations, condition 
(4.1a) implies condition (4.1b), but to be precise and mathematically rigorous, we will 
assume that both conditions hold. 

The metric for a planar gravitational system was given in the following form by 
Eardley, Liang & Sachs (1972): 

ds 2 = -d X 2 + <f> 2 (x, z) (dx 2 + dy 2 ) + V> 2 (x, z)dz 2 (4.2a) 

where 
and 

V = ?£/k{z) . (4.2c) 
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Here A = X(z), k = k(z) are independent of time. 

If A > (which is the case of physical interest — see below) , then the exact solution 
is given parametrically in terms of 77, 



X = 



A 
A 

2^ 



(sinhr/ — if) , 
(cosh, — 1) , 



A 

2^3 



A' k'\ , . / A' «'\ , (sinhr? — 7?) 

— - 2— cosh, - 1) + -— + 3— sinh, ) -i ^ 

A k / \ A « y (cosh, — 1) 



(4.3a) 

(4.36) 
(4.3c) 



where prime denotes a derivative with respect to z: ' = d/dz. It is useful to let 



(4.4) 



because at early times (small %), one recovers the standard long- wavelength initial condi- 
tions, 



= V> = X 2/3 exp j (early times) , 
which describes an initially isotropic metric. The metric variables become 

6 exp(5/) . . 



(4.5) 



X = 



<t> = 



125 (ff 

2 exp(5/) 
25 (/') 



(sinh, — ,) , 
(cosh, — 1) , 



6 exp(5/) 
125 (ff 



Bf - 2f) (cosh, - 1) + (-5f + 3 £ ) sinh, Ml| 



(4.6a) 
(4.66) 

(4.6c) 



Following Tomita (1975), we expand our results in a Taylor series in x 2 / 3 : 



, = 5 X 1/3 /' exp (-*f) 1 _ _^/3 { ff exp ( 



10/ 



25 4/3 4 / 20/ 
56 X (/) eXP V 3~ 



+ 



+ . 



(4.7) 



2/3 pYT , / ^ 



<j> = x 1 exp 



l + ^(/') 2 exp(-^/ 



(4.8) 
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ip = % 2//3 exp 



3 



(ff ) X 2/3 exp 



10/ 



+ 



(-If (/') V" + ff| (ff) X^exp 



20/ 



+ • • 



(4.9) 



It can be shown that eqs.(4.8) and (4.9) are in exact agreement with the improved spatial 
gradient expansion given in eq.(3.12). 

Let's now focus our attention on the last equation, and apply the approximations of 
eq.(4.1) which leads to: 



i> = x 



2/3 



i + ^V /3 -|(/f/V /3 + ... 



(4.10) 



When caustic formation occurs, the second term is —1 and the third term is negligible 
provided eq.(4.1b) is valid. (All higher order terms are negligible as well.) To an excellent 
approximation, the fields evolve according to 



V = 5/Y /3 

t = x 2/3 , 



V> = x 2/3 (i + lx 2/3 f") > 



(4.11) 
(4.12) 
(4.13) 



before caustic formation. The last expression is the Zel'dovich approximation because the 
corresponding density is 

4 1 

p= w\i + h 2/s f"\' (4 ' 14) 

Hence in a planar system, the Zel'dovich approximation is essentially exact until 
caustics form. Additional boundary conditions are needed to describe what happens later. 
Two possibilities are: 

(1) The particles are collisionless and they pass through each other. (These are 
cold-dark-matter particles; they comprise approximately 95% of the mass density of the 
Universe.) 

(2) The particles stick together on collision, forming a single particle whose velocity 
is given by momentum conservation. This presumably describes the baryon component 
(which comprises about 5% of the mass density of the Universe.) 

We will not discuss the issue of caustics further in this paper. 
4.2 Spherical Systems and the Improved Zel'dovich Approximation 



For spherical systems, the usual Zel'dovich approximation may not be very accurate. 
As a result, we investigate higher order corrections. 
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The Tolman-Bondi metric is 

ds 2 = -d X 2 + r 2 ( X , q) {d6 2 + sin 2 ## 2 ) + V> 2 (x, q)dq 2 , (4.15) 

where 

(|) 2 ^ + [r^)-l], (4,6a) 

* = ^/T(9), (4.166) 

and M = M(q), V = r(g) are functions of the comoving radial coordinate (/; they are 
independent of time x- (See Misner et al (1973).) 

We now choose 

M = -(/ 3 exp (5/ ((/)), 

so that for small the metric variables evolve according to 

r/q = i= = X 2/3 exp (*f) , (4.17) 



which yields an isotropic metric at early times. 

For T 2 — 1 > 0, the exact solution is analogous to the planar case 

x = 5^r 5 fp (sinh "-" ) ' (418o) 

r = T^!\f ( cosh 1 - 4) ■ (4.186) 

If T 2 — 1 < 0, then the above formulae are still valid; in this case Vr 2 - 1 and rj are pure 
imaginary. 

The planar case may be recovered if one identifies 4> = r/q, and then assumes 

q% » 1; (4.19) 
dq 

crudely speaking, this is the large radius limit. However, we have chosen to treat the 
spherically symmetric case separately because the results are qualitatively different. 
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Once again, we expand in a Taylor series in % 2 / 3 



r = X 2/ Vxp( ^) <! 1 + ^X 2/3 



X 



+ 



75 
112 

2875 
4032 



4/3 



X 



df 

dq 

df 
dq 



+ 



+ 



6_d£ 
bq dq 

6_d£ 
bq dq 



dq J bq dq 



exp 



exp 



10/ 



20/ 



exp (-10/) + 



(4.20a) 



ip = x 2 ^ 3 exp 



28 A 




2/3 



+ 



575 
224 



X 



i + x 



dq ) bq dq 
dq J bq dq 



Sd*£ 
2 dq 2 

dq 2 

(Pi 
dq 2 



5fdf 
4 \dq 

2q dq 

l_dj__ 
3q dq 



exp 



h_(dj_ 
4\dq 

25 (df 



10/ 



exp 



(4.206) 



20/ 



18 \dq 



exp (-10/) + . 



(Of course, eqs.(4.20) are in exact agreement with the improved spatial gradient expansion 
eq.(3.12).) In general, this series does not truncate after the first two terms. For example, 
consider /(g) = / cos(Ag) at q = 0. 

In Fig. 2a, for T 2 - 1< 0, we show the exact solution for the radius r(x, q) (bold line) 
with q = 0. This case describes the collapse of a high density region which is spherically 
symmetric. The light line is the first term (long- wavelength approximation). Also shown 
are the terms of order 3, 5 and 7 in spatial gradients. It is encouraging that the series 
appears to be converging before total collapse occurs. In Fig. 2b we consider T 2 — 1 > 
describing the evolution of a void region. Here, the series does not converge for x > 1.4 
(which is also the radius of convergence of Fig. 2a) . Even then the higher order terms are 
relatively close to the exact result. 



5. COMPARISON WITH THE BUCHERT-EHLERS SCHEME 



Various schemes offering higher order improvements to the Zel'dovich approxima- 
tion within Newtonian theory have been suggested, including Buchert (1989, 1992, 1993), 
Buchert & Ehlers (1993), Moutarde et al (1991), Gramann (1993), among others. Nat- 
urally one would like to be able to compare such theories with the one presented here. 
We have concentrated on the Buchert-Ehlers second order scheme (which we derived in 
section 2), because their theory is clearly and unambiguously expressed. Although there 
is a related third order scheme (Buchert 1993), the calculations required for a comparison 
are too unwieldy to discuss here. 
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All Newtonian theories are non-local because at some stage the Newton-Poisson equa- 
tion has to be solved. We therefore discuss only the restricted model described in section 
5 of their paper which depends on the scalar function f(q) which obeys the constraint 

£ijkf,jf,kmm . (^-1) 

For such a choice of f(q), the non-local Newtonian expression (2.17) becomes local: 

*\b) =Q l + I t 2/3 f^ + i* 4/3 (/,**/,* - f,kkli) ■ (5-2) 
Working in synchronous gauge their approximation implies a spatial 3-metric 

,4/3 dx (B)l 9x \b) 

=t 4/3 (Su + It^lu + t^ u ) + ft 2 / 3 / + t^T l 3 ) , (5.3) 

where 

"J~i3 = 56 (/,ij'fc/' + f,ikf,j k ~~ f' k ,(if,j)k ~ f' k ,kf,ij) ■ (5-4) 

If we set t = x an d ignore J 7 ^, then their resulting first order theory agrees perfectly with 
our improved third order theory, eq.(3.17). Both theories imply the standard Zel'dovich 
approximation expression for the density, eq.(3.19). However their second order theory 
7^ 0) does not in general agree with our improved fifth order approximation, eq.(3.18). 
Their Tij is not the same as our F^, eq.(3.18b). In particular the geodesies, which relate 
Eulerian to Lagrangian positions as a function of time, will in general differ. 

For spherical geometries, the two expressions, Tij and Fij, are indeed consistent. 
Moreover, it is straightforward to verify in general, that their traces coincide: T k k = F k k- 

(5) 

This implies that det'j^ij = det^- : at this order, both estimates of the density as 
a function of t and q will agree. The discrepancy between Tij and is still under 
investigation. 



6. DISCUSSION AND CONCLUSIONS 



Since the velocity potential plays such a prominent role in the Zel'dovich approxima- 
tion, we began our analysis with a relativistic action principle for a dust field x- I n the 
Newtonian limit, we used the non-relativistic HJ equation to solve for the velocity poten- 
tial which is expressed in terms of a Taylor series. We recovered the results of Buchert 
& Ehlers (1993). Unfortunately the final expressions for the density are nonlocal, which 
make further analytic progress difficult. 

In conflict with the principle of causality, nonlocal terms suggest that information 
may be transferred faster than the speed of light. By considering a general relativistic 
formulation, we avoid this problem. We solved Einstein's equations using a spatial gradient 
approximation. We derived the Zel'dovich approximation using the first two terms in 
our expansion. We also investigated higher order corrections. Consistent with causality, 
our expressions are local in that they involve only derivatives of the initial gravitational 
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potential. The simplicity of our final results suggest that they will be a useful tool in 
cosmology. For arbitrary initial conditions, we computed the epoch of collapse. We hope 
to apply our techniques to the interpretation of redshift surveys. 

As a check, we have compared our gradient expansion with two exact solutions of 
Einstein's equations: planar and spherical geometries. Using a Taylor series in the planar 
case, it is straightforward to show that the first two terms in the improved expansion 
are sufficient — all the higher order terms are very small. Hence in a general relativistic 
formulation, we rederive the fact that the Zel'dovich approximation is essentially exact for 
planar geometries (for a Newtonian derivation, see, e.g., Kofman (1991)). As was expected 
from the Newtonian theory, the results for spherical geometries differ quantitatively from 
the planar case. Near zero radius, the series does not truncate for a spherical system, and 
the higher corrections to the Zel'dovich approximation are indeed significant. (However, in 
the large radius limit the spherical case reduces to the planar one.) Our general relativistic 
approach, of, course is valid in general — not just for spherical and planar geometries. 

Comparison with Newtonian theories, e.g., that of Buchert and Ehlers (1993), is 
difficult because they are non-local. For the special case of spherical geometries, our 
results are completely consistent. However, in general we have not been able to show that 
two methods are equivalent. In any case, our relativistic formalism does lead to simpler 
expressions for the final results (see eq.(3.18)). 

ACKNOWLEDMENTS 

D.S.S. received partial support from the Natural Sciences and Engineering Research 
Council of Canada, and the Canadian Institute for Theoretical Astrophysics (CITA Na- 
tional Fellowship). K.M.C. was supported by the Science and Engineering Research Coun- 
cil of the UK. 

REFERENCES 

Bardeen, J.M., 1980, Phys. Rev. D, 22, 1882 
Bertschinger, E. & Jain, B., 1993, MIT preprint 

Bouchet, F.R., Juszkiewicz, R., Colombi, S. & Pellat, R., 1992, ApJ, 394, L5 

Broadhurst, T.J., Ellis, R.S., Koo, D.C. & Szalay, A.S., 1990, Nature, 343, 726 

Buchert, T., 1989, A&A, 223, 9 

Buchert, T., 1992, MNRAS, 254, 729 

Buchert, T., 1993, MNRAS, (in press) 

Buchert, T & Ehlers, J., 1993, MNRAS, 264, 375 

Comer, G.L., Deruelle, N., Langlois, D. & Parry, J., 1994, Meudon preprint 

Croudace, K.M., Parry, J., Salopek, D.S. & Stewart, J.M., 1994, ApJ, (in press) 

Dekel, A., Bertschinger, E. & Faber, S.M., 1990, ApJ, 364, 349 

De Lapparent, V., Geller, M.J., & Huchra, J.P., 1986, ApJ, 302, LI 

Ehlers, J., 1961, Abh. Math-Naturwiss. Klasse der Akad. Wiss. und Lit., Mainz, nr. 11 



20 



Eardley, D., Liang, E. & Sachs, R., 1972, J. Math. Phys., 13, 99 
Goldstein, H., 1981, Classical Mechanics, Addison Wesley 
Gramann, M., 1993, ApJ, 405, L47 

Kofman, L., 1991, in Sato, K., ed., Proc. IUPAP Conf. Primordial Nucleosynthesis and 
Evolution of the Early Universe, p. 495, Kluwer, Dordrecht 

Landau, L. & Lifshitz E.M., 1975, The Classical Theory of Fields, (fourth English edition) 
Pergamon (Oxford) 

Lifshitz, E.M. & Khalatnikov, I.M., 1964, Usp. Fiz. Nauk, 80, 391 [Sov. Phys. Usp., 6, 
495 (1964)] 

Matarrese, S., Pantano, O. & Saez, D., 1993, Phys. Rev. D, 47, 1311 

Misner, C.W., Thorne, K.S. & Wheeler, J. A., 1973, Gravitation, Freeman (New York) 

Moutarde, F., Alimi, J.-M., Bouchet, F.R., Pellat, R. k Ramani, A., 1991, ApJ, 382, 377 

Nusser, A., Dekel, A., Bertschinger, E. & Blumenthal, G.R., 1991, ApJ, 379, 6 

Parry, J., Salopek, D.S. & Stewart, J.M., 1994, Phys. Rev. D (in press) 

Peebles, P.J.E., 1980, Large Scale Structure of the Universe, (Princeton University Press) 

Peres, A., 1962, Nuovo Cim., 26, 53 

Salopek, D.S., 1992, Phys. Rev. Lett., 69, 3602 

Salopek, D.S. & Stewart, J.M., 1992, Class. Quantum Grav., 9, 1943 
Salopek, D.S. & Stewart, J.M., 1993, Phys. Rev. D, 47, 3235 

Salopek, D.S., Stewart, J.M., Croudace, K.M. & Parry, J., 1994, in Sato, K., ed., Proc. of 
37th Yamada Conference, Evolution of the Universe and its Observational Quest, Tokyo, 
Japan, (Universal Academic Press) 

Schutz, B.F., 1970, Phys. Rev. D, 2, 2762 

Schutz, B.F., 1971, Phys. Rev. D, 4, 3559 

Shandarin, S.F. & Zel'dovich, 1989, Rev. Mod. Phys., 61, 185 

Stewart, J.M., 1990, Class. Quan. Grav., 7, 1169 

Szekeres, P., 1975, Commun. Math. Phys., 41, 55 

Tomita, K., 1975, Prog. Theor. Phys., 54, 730 

Zel'dovich, Ya.B., 1970, A&A 5, 84 



21 



Figure Captions 



Figure 1. Higher order corrections to the standard Zel'dovich approximation are in 

general significant. For the improved fifth order calculation, the scale factor, a^J u = x 2 ^ 3 , 
at the time of collapse is shown as a function of the eigenvalues di of —f,ij where — / 
is the initial gravitational potential. From the bottom to the top, the curves represent 

a^coll/ a coll = 2.0,1.2,1.0 (bold), 0.9, 0.8 and 0.755 (sharp corner), where the first and last 
are the maximum and minimum values, respectively. d\ is the largest positive eigenvalue 
and a^Ju = 2/(3di) is the standard Zel'dovich approximation result. 

Figure 2. We compare our gradient expansion (1st, 3rd, 5th and 7th order) with the 
exact Tolman-Bondi solution (bold curve) for a spherical system. In (a), spherical collapse 
is shown. Apparently, the spatial gradient expansion is converging. In (b), the radius 
r expands indefinitely. The radius of convergence Xcrit = 1-4 is finite and it coincides 
with the collapse time given in (a). Even beyond the radius of convergence, the successive 
approximations do not differ much from the exact solution. 
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